Turbulent transport and heating of trace heavy ions in hot, magnetized plasmas 



M. Barnes, 1, 2, [] F. I. Parra, 1 and W. Dorland 3 

Plasma Science and Fusion Center, Massachusetts Institute of Technology, Cambridge, MA 02138, USA 
2 Oak Ridge Institute for Science and Education, Oak Ridge, TN 31831, USA 
s Department of Physics, University of Maryland, College Park, MD 20740, USA 

Scaling laws for the transport and heating of trace heavy ions in low-frequency, magnetized plasma 
turbulence are derived and compared with direct numerical simulations. The predicted dependences 
of turbulent fluxes and heating on ion charge and mass number are found to agree with numerical 
results for both stationary and differentially rotating plasmas. Heavy ion momentum transport 
is found to increase with mass, and heavy ions are found to be preferentially heated, implying a 
mass-dependent ion temperature for very weakly collisional plasmas and for partially-ionized heavy 
ions in strongly rotating plasmas. 



Introduction. Heavy ions are present in hot, magne- 
tized plasmas both in laboratory experiments and in na- 
ture. These heavy ions are often trace, i.e., their densities 
are small enough that they have only a small direct effect 
on the bulk plasma dynamics. Nonetheless, trace heavy 
ions are important in numerous contexts: main ion prop- 
erties are often inferred from heavy ion measurements be- 
cause heavy ions radiate more readily pQ; accumulation 
of heavy ions leads to dilution and increased radiative 
energy losses in magnetic confinement fusion [5J [3] ; and 
temperature measurements of minority ions in space and 
astrophysical plasmas indicate the existence of a novel 
heating mechanism [3H6]. 

Considerable effort has gone into understanding the 
particle transport of trace heavy ions, or impurities, in 
the context of magnetized toroidal plasmas for fusion. In 
particular, the scaling with charge number Z and mass 
number A of the impurity particle flux were predicted 
with a quasilinear fluid model and found to be in rela- 
tively good agreement with numerical and experimental 
results 018]. However, little to no work has been done 
on impurity momentum and energy fluxes or for turbu- 
lent heating of impurities. The latter may play a role 
not only in fusion plasmas, but also in the context of as- 
trophysical plasmas, where the temperature of minority 
ions has been observed to increase with increasing ion 
mass (4j-[6]. Cyclotron heating [9] and stochastic heating 
via large-amplitude fluctuations [TQ| have been proposed 
as possible explanations for this mass dependence. The 
turbulent heating mechanism described here provides an 
alternative explanation for the mass dependence of the 
minority ion temperature that is present even for low 
frequency, low amplitude fluctuations. 

In this Letter we use local, nonlinear, <5/-gyrokinetic 
theory [TThT3] to provide scaling predictions for trace 
heavy ion particle, momentum, and energy fluxes, as 
well as turbulent heating in hot, magnetized plasmas. 
This approach has already proven successful in determin- 
ing scalings of temperature-gradient driven turbulence in 
tokamaks |14j . We consider an inhomogeneous, axisym- 
metric plasma rotating toroidally at angular frequency 
w^, immersed in a curved, inhomogeneous magnetic field. 



To simplify our analysis, we restrict our attention to a 
region of plasma with rotation speed well below the ion 
sound speed but with a strong rotation gradient. We 
also consider only moderate values of /? — Sirp/B 2 < 1, 
where p is the mean plasma pressure and B is the mean 
magnetic field magnitude. This is directly applicable to 
toroidal confinement experiments in magnetic confine- 
ment fusion, but the scaling laws we obtain are general: 
they do not change for a stationary, homogeneous plasma 
slab and therefore also pertain to various space and as- 
trophysical plasmas. 

Gyrokinetic turbulence. The <5/-gyrokinetic theory is 
obtained by performing an asymptotic expansion in the 
small ratio of the Larmor radius, p, to system size, L, and 
averaging over the fast Larmor motion of particles. It is 
valid for low-amplitude turbulence with time scales long 
compared to the Larmor frequency, f2, and spatial scales 
comparable to p and L in the directions across and along 
the mean magnetic field, respectively. While initially de- 
veloped for magnetic confinement fusion plasmas, Sf- 
gyrokinetics can also be applied to small-scale turbu- 
lence in the solar wind, solar corona, accretion disks, and 
galaxy clusters [TH [TS] . 

We use (R, p, e) as our coordinate system, where R 
is the position of the center of a particle's Larmor or- 
bit, e = mv 2 /2 its kinetic energy, and p — mv\/2B its 
magnetic moment, with m its mass and v its speed. The 
subscripts _L and || are used to denote the components 
perpendicular and parallel to the mean magnetic field, 
respectively, with the magnetic field magnitude given by 
B. With this choice of coordinates, the electromagnetic 
gyrokinetic equation governing the evolution of the fluc- 
tuating piece of the distribution function, 5f s , is 

| + R S 'vL - (C[Sf s ]) s 
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= -( v x) s • f VF M ,a + W^^Fm.J , 

where g s = 5f s + Z s eF M>s (<& - {x) s )/T s , {.) s denotes 
an average over Larmor angle at fixed R s , (x) s ~ 
($-U||<L4||/c+ f£ 3 dp' s 5B\\/Z s e) s [26], $ is the fluctu- 
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ating electrostatic potential, SAi\ and SBu are the parallel 
components of the fluctuating magnetic vector potential 
and magnetic field, respectively, Z s is the charge num- 
ber, e the proton charge, c the speed of light, T s the 
mean temperature, Fm, s is a stationary Maxwellian dis- 
tribution of velocities in the frame rotating with velocity 
u = E?uj^V(j), <j> is the toroidal angle, R the plasma ma- 
jor radius, D/Dt = d/dt + u-S7 , R s = V| +vm,s + (v x ) s , 
with v A f. s = b/Q a x (uj?b • Vb + v 2 yB/2B) the drift 
velocity due to a curved, inhomogeneous mean magnetic 
field and (v x ) = cbx V (\) s / B the drift due to the fluc- 
tuating electromagnetic fields, b the unit vector along 
the mean magnetic field, £l s = Z s eB /m s c the Larmor 
frequency, and C describes two-particle Coulomb inter- 
actions. Plasma species is indicated by the subscript s, 
which we henceforth drop unless it is needed to avoid 
ambiguity. 

By definition, the trace ions considered here do not 
contribute to the fields. They are instead determined 
solely by the electron and main ion dynamics through 
the low-frequency Maxwell's equations, supplemented by 
the quasineutrality constraint: 

= ^Z. f d\ Sf s , (2) 
ViM|, ^-^Y^Z s eJ d 3 v v\\5f., (3) 
V ± £B|| = ^ ^ Z s e y d 3 v (b x vj.) */ a , (4) 



where $ enters Eqs. ([2|4j) through the definition for Sf 
given below Eq. (JlJ. 

With g and {$, SA\\ , 5B\\} specified by Eqs. (l]4), one 
can evaluate the turbulent heating, 

H = Ze ( X ((v,| + v M ) • V.9 - (C[Sf]))) A - H^, (5) 



and the turbulent fluxes, 



Sf(v x ) ■ Vr 



(6) 



5f(v x ) ■ Vr 



n = mR 2 {Sf(v ■ V0)(v x ) • Vr 



(7) 



— R 2 (b • V0)( Sf(y x ■ Vr)5A\\ 

C \ /A 



where r labels surfaces of constant mean pressure, (a) A = 
J d 3 r J d 3 v a/ J d 3 r is an integral over all velocity space 
and over a volume of width w (p -C w -C L) encom- 
passing the mean magnetic field line of interest, and T, 
fl, and Q are the particle, toroidal angular momentum, 
and energy fluxes, respectively. Note that the momen- 
tum flux defined in Eq. (fsl) does not include each species' 
contribution to the Maxwell stress. 



TABLE I: Scalings, S, for turbulent fluxes and heating 
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Expansion in A 1 / 2 . To obtain scaling laws for the tur- 
bulent fluxes and heating of trace heavy ions, we take 
Z ~ A I , duj^/dr ~ v t i/R 2 , and expand g = <7o+<?i+--- 
in powers of A 1 / 2 . Here Utj is the main ion thermal speed. 
We restrict our attention to (3 — %irp/B 2 < 1, and assume 
the collisional mean free path is sufficiently long that col- 
lisions may be neglected in our analysis. In what follows, 
we keep Z and A dependences separate so that we can 
consider the subsidiary expansion A x l 2 <C Z <C A. 

Because the heavy ions are trace, their space and time 
scales are those of the bulk plasma turbulence. Thus, 
Z and A only enter Eq. ([!]) through explicit factors of 
m, v ~ Vt, and Z, as well as through g itself. In what 
follows, we assume the ratio of the heavy ion to proton 
temperature is much smaller than A, giving vt ~ A~ 1/>2 . 
The two lowest order equations in our expansion are thus 
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T 



-F M (v A ) ■ i?Vw - ( V£; ) • VF M , 



where = cb x V$/£? and = uyb x V5A\\/B. 

There are two possible scalings for both go and g± 
due to a competition between terms with different A 
and Z dependences in Eqs. ^ and |To| ). In particular, 
50 oc A 1/2 , Z/A 1 / 2 , and gi oc 1, Z/A. By consider- 
(8) ing the limits |dw0/dr| < (Z/A)w ti /i? 2 and |<A^/dr| > 
(Z/A)v t i/R 2 , the number of such scalings is reduced. 
The A- and Z-scalings for g and in these limits, as 
well as for the general case, are summarized in Table [I] 

Flux and heating scalings. If go{ v \\) is a solution 
to Eq. ([9j), then —go(—v\\) is also a solution. Thus, 



dv\\ go{&, SA\\ , 5B\\ } = 0, where the overline denotes 
a statistical average. As a result, go does not contribute 
to the lowest order (i.e., electrostatic) heating or particle 
and heat fluxes, Eqs. |5|-([7|, whose integrands are oth- 
erwise even functions of v\\ ■ Conversely, the lowest order 
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FIG. 1: Normalized particle flux, (F s /n s vu)(L / pi) 2 , vs. mass 
number, A, for cases with and without differential rotation, 
cj0. The dashed line is a least-squares fit using our scaling 



predictions, given by —0.7 + 2.3/ A and 0.4- 
and right plots, respectively. 



1.6/x for the left 



FIG. 2: Normalized heat flux, {Qs/n a TiVti){L/ pi) 2 , vs. mass 
number, A, for cases with and without differential rotation, 
The dashed line is a least-squares fit using our scaling 
predictions, given by 2.0 + 15/ A and 3.7 — 3.7/^4 for the left 
and right plots, respectively. 



momentum flux integrand has a component proportional 
to mvi\, so II ~ mv t go cx A 1 / 2 go- Using our scalings for 
go, we see that II has competing terms scaling as Z and 
A, respectively. 

Note that Eq. (10) has a vu symmetry that is opposite 
that of Eq. (TqI) : if gi\V\\ ) is a solution, then gx (— vii ) is also 
a solution. Furthermore, for all higher order equations, 
one can show that the symmetry in v\\ alternates between 
that of Eqs. ([9]) and (10). As a result, the only compo- 
nents of g that contribute to the particle and heat fluxes 
and heating are g\, #3, etc. Using Eqs. ^ and ([7|, we 
have {T,Q} ~ g\, which in the general case has compet- 
ing terms scaling as ZJA and 1 (no Z or A dependence), 
respectively. 

The first term in the heating expression, ([5]), is the 
Joule heating and is scaled up by an explicit factor of 
Z (arising from the current), while the second term is 
viscous heating. At lowest order, the Joule heating term 
gives H cx {Zv\\g Q , Zgi), giving H cx (Z, Z 2 /A). The 
viscous heating is proportional to II cx (Z, A). H thus 
has competing terms scaling as Z 2 /A, Z, and A, respec- 
tively. The scalings of the various fluxes and heating are 
summarized in Table |U 

Minority ion temperature. Integrating Eq. ^ by 
parts in time and using Eq. ([!]), the heating can be ex- 
pressed as [H ED EE] 



T s Sf s 

Fm,s 



C[Sf. 
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dlnT s 
dr 



<91nn s 
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(11) 



Our scalings indicate that H s increases in magnitude 
with A or Z, but T s and Q s do not. The first term 
in Eq. (11) must thus dominate for A or Z large. This 



term, which we identify as the collisional entropy gener- 



ation, is positive definite when summed over species. We 
argue that it is also positive species by species for the low 
collisionalities considered here. 

The collision operator, C, consists of a test-particle 
piece, which is a diffusion operator in velocity space, and 
a field-particle piece, which is an integral operator |19) . 
Both contributions are inversely proportional to the col- 
lisional mean free path and thus small, except at small 
scales in the velocity space where large derivatives in the 
test-particle operator compensate [TBI HOI 121] ■ The test- 
particle operator should thus dominate in weakly colli- 
sional plasmas, and its diffusive nature ensures that its 
contribution to entropy generation is positive-definite. 

Consequently, trace heavy ions must be heated by tur- 
bulence instead of cooled. For this heating process to 
subside, the trace ion temperature must become large 
enough to interfere with our large A expansion. This 
happens when the heavy ion temperature exceeds the 
main ion temperature by a factor of A ~ Z. In this 
limit, the turbulent heating H becomes comparable to 
the heat flux Q so that H is no longer required to be 
positive definite. Our theory thus predicts that heavy 
ions will be hotter than light ions by a factor of A ~ Z 
- but only if turbulent heating is larger than collisional 
temperature equilibration. 

The collisional temperature equilibration of the main 
ions, i, and a trace heavy ion species, s, is £ s = 
(8/3^)(Z 2 /A s )n s AT s v tl /X mfp , where AT S = T. - T u 
and A m f p is the mean free path for collisions between 
the main ions. From Eq. ([5]), we estimate H s ~ 
Sn s Ti(Srii/ni) 2 Vti/ L, where S is the scaling of H with 
A and Z given in Table |lj and we have assumed e<&/Ti ~ 
Srii/rii ~ 5n s /n s . The ratio of turbulent heating to 
collisional temperature equilibration is thus H/S ~ 
5(^/Z 2 )(A mfp /L)(fe l /n l ) 2 ~ S{A/Z 2 )\ miv /v tl T El with 
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FIG. 3: Normalized toroidal angular momentum flux, 
(J\ s / m,in s Lvti){L / pi) 2 , vs. charge number, Z, and mass num- 
ber, A, for cases with and without differential rotation, cj^. 
The dashed lines are least-squares fits using our scaling pre- 
dictions, given by and —1.7 A — 5.1 for the left and right 
plots, respectively. The fact that II — for the case with 
duj^/dr = is a consequence of a symmetry property of the 
gyrokinetic equation [25] . 



te the characteristic time scale over which the equilib- 
rium density and temperature vary. 

Numerical results. To test our predictions for the 
scalings of turbulent transport and heating, we em- 
ploy the local, 8 /-gyrokinetic code GS2 [25] . We con- 
sider an axisymmetric system with sheared magnetic field 
lines mapping out nested toroidal surfaces with circu- 
lar cross sections (known as the Cyclone Base Case [23] 
and parametrized using the Miller local equilibrium 
model [H]). Each simulation is electrostatic and includes 
kinetic electrons, as well as kinetic main and trace heavy 
ions with a wide range of Z and A values. The turbu- 
lence is driven by gradients in the mean ion and electron 
densities and temperatures, with i? (dlnn/ 'dr) = 2.2 for 
the electrons and main ions, and Ro(dlnT/dr) — 6.9 for 
all species, with R the major radius at the center of the 
constant pressure surface. The collision frequency is cho- 
sen small, Ro/X m { p = 0.003, so that heavy ion collisions 
do not affect our scalings. 

Two sets of simulations were carried out: one with a 
stationary plasma (dui^ /dr — 0) and one with a differen- 
tially rotating plasma {duj^/dr — 4.67v t i/ R 2 ). The simu- 
lation results are shown in Figs. Q-Q. Data points for 
fluxes and heating at various Z and A values are plotted 
as solid circles and fit using a least-squares analysis with 
the predicted lowest order Z and A dependences, as well 
as the first order correction. In each case, the predicted 
scalings fit the data well. It should be noted that the 
momentum flux for dw^, /dr = is zero for all species due 
to a fundamental symmetry of the (^/-gyrokinetic equa- 
tion [25] . 

Discussion. We now discuss the implications of the 
trace heavy ion scalings derived in this Letter. First, the 



FIG. 4: Normalized heating, H s (L/n s TiV t i){L/ p^ 2 , 
vs. charge number, Z, and mass number, A, for cases 
with and without differential rotation, The dashed lines 
are least-squares fits using our scaling predictions, given by 
1.1Z 2 /A - 1.8 and 1.1A + 5.0 for the left and right plots, 
respectively. 



preferential heating of heavy ions should lead to large 
temperature disparities between different ion species in 
nearly collisionless plasmas. Many space and astrophysi- 
cal plasmas are weakly collisional enough (i.e., (8n/n) 2 > 
L/\ m f p ) that turbulent heating should dominate over col- 
lisional equilibration, and preferential heating of heavy 
ions is indeed observed |U [S]. However, for such low 
collisionalities the equilibrium can deviate strongly from 
the isotropic Maxwellian assumed in our analysis, which 
cannot consequently address the large 2j_/Tm values ob- 
served in coronal holes and the fast solar wind [5J. 

Magnetic confinement fusion plasmas typically have 
(8n/n) 2 < L/X m [ p so that collisional temperature equi- 
libration dominates over turbulent heating and all ions 
have the same temperature. However, for rotating plas- 
mas our results indicate that the turbulent heating is 
enhanced by an additional factor of (A/Z) 2 relative to 
the equilibration. It may therefore be possible for heavy, 
partially ionized impurities to be heated by turbulence 
to temperatures significantly larger than the main ions. 

Because the momentum transport of heavy ions is en- 
hanced by A, it can generate flows of order the ion ther- 
mal speed for densities as small as rii/A. In this limit, 
heavy ions could thus significantly alter bulk plasma mo- 
mentum transport. 

We thank S. C. Cowley, E. Quataert, and A. A. 
Schekochihin for useful discussions. M.B. was supported 
by a US DoE FES Postdoctoral Fellowship, F.I.P. was 
supported by US DoE Grant No DE-FG02-91ER-54109, 
and computing time was provided by HPC-FF (Jiilich). 



Electronic address: mabarnes@mit.edu 



5 



[1] R. C. Isler, Plasma Phys. Control. Fusion 36, 171 (1994). 
[2] E. B. Meservey, N. Bretz, D. L. Dimock, and E. Hinnov, 

Nucl. Fusion 16, 593 (1976). 
[3] M. Z. Tokar el al, Nucl. Fusion 37, 1691 (1997). 
[4] W. K. H. Schmidt, H. Rosenbauer, E. G. Shelly, and J. 

Geiss, Geophys. Res. Lett. 7, 697 (1980). 
[5] M. R. Collier et al, Geophys. Res. Lett. 23, 1191 (1996). 
[6] J. L. Kohl et al, Solar Physics 175, 613 (1997). 
[7] C. Angioni and A. G. Peeters, Phys. Rev. Lett. 96, 

095003 (2006). 

[8] C. Angioni et al, Plasma Phys. Control. Fusion 51, 
124017 (2009). 

[9] S. R. Cranmer, G. B. Field, and J. L. Kohl, Astrophys. 

J. 518, 937 (1999). 
[10] B. D. G. Chandran et al, Astrophys. J. 720, 503 (2010). 
[11] P. J. Catto, Plasma Phys. 20, 719 (1978). 
[12] E. A. Frieman and L. Chen, Phys. Fluids 25, 502 (1982). 
[13] G. G. Howes et al, Astrophys. J. 651, 590 (2006). 
[14] M. Barnes, F. I. Parra, and A. A. Schekochihin, Phys. 

Rev. Lett. 107, 115003 (2011), |arXiv:1104.4514| 
[15] G. G. Howes et al, J. Geophys. Res. 113, A05103 (2008). 
[16] A. A. Schekochihin et al, Astrophys. J. Suppl 182, 310 

(2009), arXiv: 0704.0044. 



[17] J. A. Krommes and G. Hu, Phys. Plasmas 1, 3211 (1994). 
[18] I. G. Abel et al, Reviews on Progress in Physics, sub- 
mitted (2012). 

[19] P. Helander and D. J. Sigmar, Collisional transport in 
magnetized plasmas (Cambridge University Press, 2002). 

[20] I. G. Abel et al, "Phys. Plasmas 15, 122509 (2008), 
larXiv:0806. 10691 

[21] M. Barnes et al, Phys. Plasmas 16, 072107 (2009). 

[22] W. Dorland, F. Jenko, M. Kotschenreuther, and B. N. 
Rogers, Phys. Rev. Lett 85, 5579 (2000). 

[23] A. M. Dimits et al, Phys. Plasmas 7, 969 (2000). 

[24] R. L. Miller et al, Phys. Plasmas 5, 973 (1998). 

[25] F. I. Parra, M. Barnes, and A. G. Peeters, Phys. Plasmas 
18, 062501 (2011). 

[26] The fields <£, 6 An, and <5B|| are independent of Lar- 
mor angle at fixed particle position, r, but not at fixed 
R = r + vi x b/fi. Thus care must be taken to specify 
which spatial coordinate is held fixed for velocity inte- 
gration. The /i-integral contained in {\) is performed at 
fixed R, but all other velocity integrals in this Letter are 
performed at fixed r. 



